Fluctuating Glasma initial conditions and flow in heavy ion collisions 
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We compute initial conditions in heavy-ion collisions within the Color Glass Condensate (CGC) 
framework by combining the impact parameter dependent saturation model (IP-Sat) with the clas- 
sical Yang-Mills description of initial Glasma fields. In addition to fluctuations of nucleon positions, 
this IP-Glasma description includes quantum fluctuations of color charges on the length-scale de- 
termined by the inverse nuclear saturation scale Qs. The model naturally produces initial energy 
fluctuations that are described by a negative binomial distribution. The ratio of triangularity to 
eccentricity £3/62 is close to that in a model tuned to reproduce experimental flow data. We com- 
pare transverse momentum spectra and U2,3,4(pt) of pious from different models of initial conditions 
using relativistic viscous hydrodynamic evolution. 



A large uncertainty in the hydrodynamical description 
of ultrarelativistic heavy ion collisions is our imperfect 
knowledge of multi-gluon states in the nuclear wavefunc- 
tions and the early-time dynamics of gluon fields after the 
collision. Recent studies of observables sensitive to har- 
monics of hydrodynamic flow distributions provide con- 
straints both on the shear viscosity of the Quark-Gluon 
Plasma (QGP) and the initial state dynamics [1-12]. In 
particular, higher flow harmonics are sensitive to quan- 
tum fluctuations in multi-particle production and dy- 
namics. 

An ab initio framework for multi-particle production 
is the Color Glass Condensate (CGC) [13] wherein the 
initial state dynamics is described by flowing Glasma 
fields [14, 15]. There are several sources of quantum 
fiuctuations that can infiuence hydrodynamic flow on an 
event- by- event basis. An important source of fluctua- 
tions, generic to all models of quantum fluctuations, arc 
fluctuations in the distributions of nucleons in the nu- 
clear wavefunctions. In addition there are fluctuations 
in the color charge distributions inside a nucleon. This, 
combined with Lorentz contraction, results in "lumpy" 
transverse projections of color charge configurations that 
vary event to event ^. The scale of this lumpiness is given 
on average by the nuclear saturation scale Qs, which 
corresponds to distance scales smaller than the nucleon 
size [18] . For each such configuration of color charges, the 
QCD parton model predicts dynamical event-by-event 
fiuctuations in the multiplicities, the impact parameters 
and the rapidities of produced gluons [19]. 

All these sources of fluctuations are captured in the 
CGC Glasma flux tube picture^ of multi-particle pro- 
duction [21-23]. The relevant feature of this scenario 



^ Further "pomeron loop" fluctuations arc suppressed in large nu- 
clei [16] though they may be signiflcant in nucleon-nucleon colli- 
sions [17]. 

^ Theoretical support for this framework is obtained from factor- 
ization theorems for inclusive multi-particle production in A+A 
collisions [20]. 



is that long range rapidity correlations from the initial 
state wavefunctions, coherent over 1/Qs in the transverse 
plane, are efficiently transmitted into hydrodynamic flow 
of the final state quark-gluon matter [24, 25]. 

Recently, JVIonte-Carlo Glaubcr-typc models (JVIC- 
Glauber) and IVIonte-Carlo implementations of the KLN 
model (JMC-KLN) [26, 27] have been compared to ex- 
perimental data on elliptic and triangular moments of 
the fiow distribution. While both types of models treat 
fiuctuations in nucleon positions identically, the Glauber 
model implementations do not specify a mechanism for 
multi-particle production which would constrain the ini- 
tial energy density distribution. JMC-Glauber initial con- 
ditions [7, 28] can be tuned to reproduce data on both 
elliptic and triangular flow from RHIC and the LHC. 
The MC-KLN model is motivated by the CGC with ap- 
proximations that will be discussed further below. It re- 
quires larger values of the viscosity to entropy density ra- 
tio (ry/s) relative to the Glauber model values to describe 
elliptic flow data. This however leads it to underpredict 
triangular flow data. 

Odd flow harmonics are entirely driven by fluctua- 
tions; it is therefore essential to have a realistic descrip- 
tion of quantum fluctuations in multi-particle produc- 
tion to properly describe the final state dynamics. To- 
wards this end, we will consider in this letter the impact 
parameter dependent saturation model (IP-Sat) [29, 30] 
to determine fluctuating configurations of color charges 
in two incoming highly energetic nuclei. This model 
is formally similar to the classical CGC JVIcLerran- 
Venugopalan (JVIV) model of nuclear wavefunctions [31], 
but additionally includes Bjorken x and impact parame- 
ter dependence through eikonalized gluon distributions of 
the proton that are constrained'^ by HERA inclusive and 
diffractive e-l-p deeply inelastic scattering (DIS) data [32]. 



^ The IP-Sat model gives good x-squared fits to available small x 
HERA data [32] and fixed target e-|-A DIS data [18]. Since the 
analysis of Ref. [32] , more precise data is now available; a repeat 
analysis is desirable. 
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Most importantly, the model is in excellent agreement 
with data on n-particle multiplicity distributions in p+p 
collisions at RHIC and the LHC and in A+A collisions at 
RHIC [33, 34], an essential requirement for microscopic 
models. The MC-KLN model does not contain these fea- 
tures; a scheme to introduce fluctuations in the model 
has only been discussed recently [35]. 

The color charges p''{x~,xi_) in the IP-Sat/MV 
model act as local sources for small x classical gluon 
fields. These are determined by solving the CYM equa- 
tions [Di^i,F'^'^] = , with the color current J^^ = 
(5^+P(^)(x~, xx) generated by a nucleus A (B) moving 
along the a;"*" {x~) direction. The solution in light cone 
gauge A^{A~) = are the pure gauge fields [31, 36, 37] 
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path ordered Wilson line in the fundamental represen- 
tation, where the infrared regulator m (of order Aqcd) 
incorporates color confinement at the nucleon level. The 
solution of the CYM equations in Schwinger gauge = 
for the two nuclei at r = has a simple expression in 
terms of the gauge fields of the colliding nuclei [14]: 
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(1) 



and'' drA' = 0, a^A'' = 0. In the limit r 0, 
A^ = —Eri/2, with the longitudinal component of 
the electric field. At t = 0, the only non-zero compo- 
nents of the field strength tensor are the longitudinal 
magnetic and electric fields, which can be computed non- 
perturbatively. They determine the energy density of the 
Glasma at r = at each transverse position in a single 
event [15, 38]. The Glasma distribution computed from 
the CYM equations^ (IP-Glasma) is matched event-by- 
event to viscous relativistic hydrodynamics [4, 5] to com- 
pute harmonics of the fiow distributions. 

We will now discuss details of the computation. Nu- 
cleon positions in the nucleus are sampled from a Fermi 
distribution. The saturation scale bj^) is de- 

termined from the IP- Sat dipole cross section for each 
nucleon, where bj_ is the impact parameter relative 
to each nucleon's center. The color charge squared 
per unit transverse area (7^/i^(x,bx) is proportional^ 
to ^p^(x,b_L). For the nucleus, g'^ ^i\{x,-x.±_) is ob- 
tained [41] by adding the individual nucleon g^/x^ at the 
same x and transverse position x_l in the nucleus. 



The metric 
diag(l, -1, - 



the (r, xx,»?) coordinate system is = 
-T-^) so that Arj = —t'^A^. The ± components 



of the gauge field are related by = itx^A''. 
As noted previously [39], the CYM approach treats soft modes 
with fcx < Qs more accurately than in commonly used k± fac- 
torized descriptions. 

The exact numerical factor between the two quantities depends 
on the details of the calculation [40] but will not be relevant for 
our final results. 
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FIG. 1. The IP-Glasma event-by-event distribution in energy 
for 6 = 9 fm on the the lattice compared to different functional 
forms. The negative binomial distribution (NBD) gives the 
best fit. 

The lattice formulation of the Glasma initial conditions 
in Eq. (1) was first given in [42]. On a transverse lattice, 
random color charges^ p°'{x±) are sampled from 



(p?(x^)p?(yx))=<5"V^'5^(x^-y^) 



g>i(x±) 



(2) 



where the indices® k, I 



1,2,..., Ny represent a dis- 
cretized x~ coordinate [40]. The path ordered Wilson 
line is discretized as 



Va(b)(x±) = n cxp -ii 
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To each lattice site j we assign two SU(A^c) matri- 



ces V( 



(A) J 



and V(B).j, each of which defines a pure 



gauge configuration with the link variables U, 



{A,B),o 



one lattice site in the i 



, where +ei indicates a shift from j by 
1, 2 transverse direction. The 
link variables in the future lightcone C/j, are determined 
from solutions of the lattice CYM equations at r = 0, 



tr^" [{ulA) + UlB))a + U''^) 



(l + C/^)(C/(% + C/, 



= 0. 



(4) 



where arc the generators of SU{Nc) in the fundamental 
representation (The cell index j is omitted here). The 
— 1 equations (4) arc highly non-linear and for N,. = 3 
arc solved itcratively. 



^ Here, and henceforth, the distributions are evaluated at a; = 
(p± ) / for zero rapidity, where {p± ) is the average transverse 
momentum of charged hadrons in p+p collisions at a given y/s. 

* Ny = 100 in all computations presented here. 
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The total energy density on the lattice at r = is given 

by 

£{T = 0) = ^iN,-RetTUa) + ^tTE^, (5) 

where the first term is the longitudinal magnetic energy, 
with the plaquette given by t/i — Uf U^,^ U^l,- U^'^ . 
The explicit lattice expression for the longitudinal elec- 
tric field in the second term can be found in Refs. [42, 43] . 
In Fig. (1) we show the evcnt-by-event fluctuation in 
the energy per unit rapidity at time r = 0.4 fm. The 
mean was adjusted to reproduce particle multiplicities 
after hydrodynamic evolution. This and all following re- 
sults are for Au+Au collisions at RHIC energies {y/s = 
200 A GeV) at midrapidity. The best fit is given by a neg- 
ative binomial (NBD) distribution, as predicted in the 
Glasma flux tube framework [44] ; our result adds further 
confirmation to a previous non-perturbative study [23]. 
The fact that the Glasma NBD distribution fits p+p 
multiplicity distributions over RHIC and LHC ener- 
gies [33, 34] lends confidence that our picture includes 
fluctuations properly. 

We now show the energy density distribution in the 
transverse plane in Fig. (2). We compare to the MC-KLN 
model and to an MC-Glauber model that was tuned to 
reproduce experimental data [4, 11]. In the latter, for 
every participant nucleon, a Gaussian distributed energy 
density is added. Its parameters are the same for ev- 
ery nucleon in every event, with the width chosen to be 
0.4 fm to best describe anisotropic flow data. We will 
also present results for a model where the same Gaus- 
sians are assigned to each binary collision. The resulting 
initial energy densities differ significantly. In particu- 
lar, fluctuations in the present computation occur on the 
length-scale Qj'^{x±), leading to finer structures in the 
initial energy density relative to the other models. As 
noted in [35], this feature of CGC physics is missing in 
the MC-KLN model. 

We next determine the participant ellipticity £2 and 
triangularity £3 of all models. Final fiow of hadrons is 
to good approximation proportional to the respective £„ 
[45] , which makes these eccentricities a good indicator of 
what to expect for z;„. We define 

■v/ (r" cos(Ti(j))y^ + (r" sin(n0))^ 
s« = p) ' (6) 

where (■) is the energy density weighted average. The re- 
sults from averages over ^ 600 events for each point plot- 
ted are shown in Fig. 3. The ellipticity is largest in the 
MC-KLN model and smallest in the MC-Glauber model 
with participant scaling of the energy density (iVpart)- 
The result of the present calculation lies in between, 
agreeing surprisingly well with the MC-Glauber model 
using binary collision scaling (A'binary)- This confirms 
previous results in the GYM framework using average 
initial conditions [46]. 




FIG. 2. (Color online) Initial energy density (arbitrary units) 
in the transverse plane in three different heavy-ion collision 
events: from top to bottom, IP-Glasma, MC-KLN and MC- 
Glauber [11] models. 

The triangularities are very similar, with the MC-KLN 
result being below the other models for most impact pa- 
rameters. Again, the present calculation is closest to the 
MC-Glauber model with binary collision scaling. There 
is no parameter dependence of eccentricities and trian- 
gularities in the IP-Glasma results shown in Fig. 3. It 
is reassuring that both are close to those from the MC- 
Glauber model because the latter is tuned to reproduce 
data even though it does not have dynamical QCD fluc- 
tuations. 

We have checked that our results for £2, £3 are insensi- 
tive to the choice of the lattice spacing a, despite a log- 
arithmic ultraviolet divergence of the energy density at 
T = [47] . They are furthermore insensitive to the choice 
of g, the ratio g^^/Qs, and the uncertainty in Bjorken x 
at a given energy. 

Finally, in Fig. 4 we present results for the transverse 
momentum spectrum and anisotropic flow of thermal 
pions after evolution using MUSIC [4, 48] with boost- 
invariant initial conditions and shear viscosity to entropy 
density ratio 77/s = 0.08. Average maximal energy densi- 
ties of all models were normalized to assure similar flnal 
multiplicities. More pronounced hot spots lead to harder 
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FIG. 3. (Color online) Average participant ellipticity (upper 
panel) and triangularity (lower panel) of the initial state. This 
calculation (circles), MC-KLN (squares), Glauber implemen- 
tation with participant and binary collision scaling (triangles). 
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FIG. 4. (Color online) Thermal 7r+ transverse momentum 
spectra (upper) and anisotropic flow coefiicients vi, va, and 
ii4 as functions of pr (lower) from IP-Glasma initial conditions 
(solid), MC-KLN (dashed), MC-Glauber using participant 
scaling (dotted) and binary collision scaling (dash-dotted). 



momentum spectra in the present calculation compared 
to MC-KLN and MC-Glauber models. Differences in 
V2{pt) a-iid V3{pt) are as expected from the initial ec- 
centricities of the different models. 

As discussed at the outset, MC-KLN fails to describe 
experimental V2 and W3 simultaneously [7, 28] because of 
its small ratio 63/62- The fluctuating IP-Glasma initial 
state presented here has a larger 83/62, closer to that of 
the MC-Glauber model that is tuned to describe experi- 
mental Vn reasonably well [11]. 

In summary, we introduced the IP-Glasma model 
of fluctuating initial conditions for heavy-ion collisions. 
This model goes beyond the MC-KLN implementation 
by using CYM solutions instead of fc^-factorization and 
including quantum fluctuations on the dynamically gen- 
erated transverse length scale l/Qs- Further, unlike MC- 
KLN, its parameters are fixed by HERA inclusive and 
diffractive e-|-p DIS data. At fixed impact parameter, this 
model naturally produces NBD multiplicity fluctuations 
that are known to describe p + p and A + A multiplicity 
distributions, and its ratio of initial triangularity to ec- 
centricity is more compatible with experimental data of 
harmonic flow coefficients. 

Looking forward, an improved matching to the hydro- 
dynamic description, starting at time tq, can be achieved 



by including classical Yang-Mills evolution of the system 
up to this time. However, we do not expect a signifi- 
cant modification of the presented results for £2 and £3 
as suggested by previous work [46]. Further refinements 
include treating color charge correlations encoded in the 
JIMWLK hierarchy for improved rapidity and energy dis- 
tributions [49, 50] and eliminating arbitrariness in choice 
of thermalization time by an ab initio treatment of ther- 
malization [51-54]. 
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